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Abstract. We present a filtered backprojection algorithm for reconstructing the 
Wigner function of a system of large angular momentum j from Stern-Gerlach- 
type measurements. Our method is advantageous over the full determination of 
the density matrix in that it is insensitive to experimental fluctuations in j, and 
allows for a natural elimination of high-frequency noise in the Wigner function 
by taking into account the experimental uncertainties in the determination of j, 
its projection m, and the quantization axis orientation. No data binning and 
no arbitrary smoothing parameters are necessary in this reconstruction. Using 
recently published data [Riedel et ah, Nature 464:1170 (2010)] we reconstruct the 
Wigner function of a spin-squeezed state of a Bose-Einstein condensate of about 
1250 atoms, demonstrating that measurements along quantization axes lying in 
a single plane are sufficient for performing this tomographic reconstruction. Our 
method does not guarantee positivity of the reconstructed density matrix in the 
presence of experimental noise, which is a general limitation of backprojection 
algorithms. 



1. Introduction 

The reconstruction of the quantum-mechanical state of a system from measurements 
is an important topic of the emerging field of quantum technology pQ . Through partial 
or full state reconstruction we can estimate entanglement properties of multipartite 
quantum systems, and judge their usefulness for further experimental progress in 
fields such as quantum metrology [21 [3l [4j [5l [6l [Tj [Sj [9] , quantum simulation 10J , and 
quantum computation [TT | \12 \ [T3l [Ml ITS] . 

Particularly in quantum metrology, experiments often involve large numbers 
of particles, and single-particle resolution is unavailable in both control and 
measurement. Because of this limitation, standard methods for reconstructing the 
quantum-mechanical density matrix |16l I17L 113] cannot be applied. For instance, and 
centrally to this work, in a Bose-Einstein condensate consisting of ./V atoms, with each 
atom representing a pseudo-spin-1/2 subsystem, the total spin length j = N/2 can take 
on very large values and the known reconstruction procedures become problematic. 
In a single Stern-Gerlach measurement on the atomic ensemble we measure the 
numbers of up and down spins N^- and iVj., in terms of which the total spin is 
j = (iV-f + ATj.) /2 and the projection quantum number is m = (iV-j- — N\)/2. Since it is 
very difficult to determine the populations N-f and N± with atomic accuracy [HI [19] , 
the density matrix, which requires knowledge of j, becomes impossible to reconstruct 
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in full. Further, reconstructing the (2j + l) 2 degrees of freedom of the density 
matrix [TnittZHHI] requires at least as many uncorrelated measurements, and therefore 
the experimental uncertainty in m will hinder this full determination. In the absence 
of reliable data, there will be significant uncertainty and noise throughout the density 
matrix in its Dicke representation p mm i — (jm\p\jm'}, which severely limits its 
usefulness. We need a method for calculating those components of p which are 
significant even in the presence of noise and for very large values of j, and a way 
of determining which components must remain unknown. 

The Wigner function [3T] is ideal for such a controlled reconstruction. It is a real- 
valued function on a sphere of radius fry/ j(j + 1), represented in terms of orthonormal 
Laplace spherical harmonics as [22] 

2j fe 

W{0,<p) = J2 J2 PkgW*>¥>)> (!) 

k=0 q=—k 

where i? is the polar angle measured from the +z axis, and tp is the azimuthal 
angle around the z axis. While this sphere is commonly called a generalized Bloch 
sphere [4 , its surface actually represents a two-dimensional phase space instead of 
a Hilbert space as for the original Bloch sphere. This Wigner function contains the 
same information as the density matrix for any spin-j system. While the marginals 
of the better-known Wigner function in planar space [21] [231 12H HI] are real-space 
or momentum-space probability distributions, the marginals of the spherical Wigner 
function are the projection quantum number distributions along all quantization axes 
[see ([6| below]; further, the expectation value of the angular momentum vector is 
proportional to the "center of mass" of the Wigner function, {(S x ), (S y ), (S z )} = 

y ■ 7 ^ +1 ^ 2j+1 ^ x J* sin i?di9 d^jsin?? cosp, sin^sintp, cos-djW^, tp). 

Most importantly, the Wigner function allows us to differentiate between more 
significant components pu q (with smaller values of k) and more noise-prone components 
(with larger values of k) in a natural way. Further, if only components with k <C 2j 
are reconstructed, then accurate knowledge of j is not necessary. As detailed in 
section [2j the transformation from j-space (the Dicke representation p mm < of the 
density matrix) to fc-space (the spherical harmonic decomposition pk q of the Wigner 
function) proceeds though coupling coefficients which, at low k, are smooth in both 
j and m; this significantly reduces the impact of uncertainties in the experimental 
determination of (j, m). 

Methods for reconstructing planar Wigner functions by inverse Radon transform 
are well established in the context of nonlinear optics (24] [25]. In the past they 
have also been applied to tomographic data on large-spin quantum systems, locally 
approximating the Bloch sphere by a tangental plane and neglecting its curvature [5] . 
While this approximation is valid for spin states which are very localized on the 
Bloch sphere and do not wrap around it, future experimental progress is expected to 
produce ever more delocalized states {e.g., Schrodinger-cat states) whose properties 
are strongly influenced by the spherical shape of the Bloch sphere. Previous work on 
the reconstruction of the Wigner function on the full Bloch sphere has used the Husimi- 
Q distribution as input [26 , which is the convolution of the system's Wigner function 
with that of a coherent state (see section [3]). This convolution washes out features 
of the Wigner function that are smaller than a coherent state. Since the principal 
characteristic of spin-squeezed states is that their Wigner function possesses a peak 
width smaller than that of a coherent state, such a deconvolution-based reconstruction 



Tomographic reconstruction of the Wigner function on the Bloch sphere 



3 



approach is ill suited for studying spin-squeezed states, which is the goal of much 
current research in atomic physics [U[31[3J|S1IH1[I1[E1[S]- ^ e therefore require a 
new method for reconstructing the complex quantum-mechanical states of large-spin 
systems from experimental data in the absence of simplifying circumstances, such as 
strong phase-space localization and/or lack of spin squeezing. 

This paper is organized as follows. In section [2] we present a novel filtered 
backprojection algorithm for reconstructing the Wigner function from experimental 
Stern-Gerlach data. Section [3] specializes this algorithm to data acquired with 
quantization axes lying in a single plane. Finally, section[4]applies the latter algorithm 
to a data set acquired in our group [5J. In what follows, a "single Stern-Gerlach 
measurement" describes a single determination of the projection quantum number m 
of a quantum system along a certain quantization axis. In our case this corresponds 
to a single run of state preparation and population determination of a two-component 
Bose-Einstein condensate, yielding a single tuple (j„,m n ). The equivalent for the 
original experiment |27) is sending a single silver atom through the experimental 
apparatus, and determining its deflection by the magnetic field gradient. On the 
other hand, a "Stern-Gerlach experiment" is a series of many single Stern-Gerlach 
measurements with fixed quantization axis, sufficient to determine the probability 
distribution {p~j,P-j+i, ■ ■ ■ ,Pj} while j is presumed fixed. 



2. Wigner function reconstruction by filtered backprojection 

The density matrix p of a system of total angular momentum j (assumed fixed here; 



this condition will be relaxed in section 2.1) is usually expressed in one of the two 
forms 

2j k 

P mm > = (jm\p\jm') = ^2^2 P^k™" 1 ( 2a ) 

fc=0 q=-k 

3 j 

P k 1 = Prnm't 3 ^" 1 , (26) 

m=-j m'=-j 

with the transformation coefficients (in the following simply termed Clebsch-Gordan 
coefficients) [55] 

47 m ' = (-iy- m -«{j,m;j,-m%q}, (3) 
nonzero only if q = m — ml '. Both forms contain the same information and are 



completely interchangeable. While form (2a) is more common, form (26) allows 
expressing the Wigner function on the Bloch sphere 0. Since our goal is the 



reconstruction of the Wigner function from experimental data, we focus on form (26), 
in particular its low-fc components. 

In order to determine the unknown quantum-mechanical state of a system of 
total spin j, it is necessary that many instances of this state can be generated 
experimentally [I], on which destructive measurements are performed. Further, 
projective Stern-Gerlach measurements must be performed along many different 
quantization axis orientations (i9„,<p„). For the correctness of the following 
reconstruction method it is crucial that these measured quantization axes are 
distributed as evenly as possible over the hemisphere of orientations. Since this 
requirement may be difficult to fulfill experimentally, we assign weights c„ to 
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the individual measurements in order for the weighted measurement density to 
approximate a homogeneous distribution of quantization axes as best possible. 
Notice that these weights are independent of the outcomes m n of the Stern-Gerlach 
measurements. In the ideal case of homogeneously distributed quantization axis 
orientations (for example through the vertices of a geodesic hemisphere), all these 
weights are chosen equal and the data are used most efficiently. 

In this way, the results from M single Stern-Gerlach measurements along various 
quantization axes orientations are assembled into a data set of tuples (p n ,c n , m n ) 
with n = 1...M and S n =i c » = Our filtered backprojection algorithm for 
reconstructing the Wigner function coefficients is then given by 

M 

= (2k + l)J2cnD k q0 (v n ,KMo nm "i (4) 

71=1 

with D 3 , {a, /3, 7) = (jm'\e~ laJz e~ ll3Jy e~~ nJ * \jm) a Wigner rotation matrix [25] ; 
in particular D q0 (tp, 0) = <J sf^i^kg^' f)- This is formally equivalent to the 
filtered backprojection algorithm used for planar inverse Radon transforms [29] . 
with the factor 2k + 1 representing the "filter", and the summand representing the 
backprojection. Our algorithm has all of the typical properties of planar inverse Radon 
transforms by filtered backprojection; no data binning is required, and there are no ad 
hoc parameters to be chosen or optimized. Further, as the backprojection algorithm 
is a direct sum and does not include an inversion step (such as a straight inversion 
of the Radon transform would require) , the impact of experimental noise is bounded 
in the result. It is this last property which makes backprojection algorithms fast and 
reliable in practical applications such as X-ray computed tomography [29] . 

Our specific backprojection Q can be interpreted in an intuitive way. The 
measured values of m n in the coordinate frame attached to the quantization axis 
{d n ,ip n ) are distributed according to the diagonal elements /? m m „ and are converted 

L kq' 



from j-space into /c-space via the Clebsch-Gordan coefficients tj™i nmn with q' = (see 



section Appendix A for a numerical procedure). They are then rotated into the lab 
frame through the rotation matrices D qq ,((p n , , d n ,Xn) with the value of Xn irrelevant 
(set to zero) since q' = 0. 

In the following, we demonstrate that this algorithm Q works in the limit of 
infinite data. If all quantization axis orientations have been used with equal frequency, 
and infinitely many measurements have been performed along each quantization 
axis, the sum over measurements X^^Li c " can ^ e replaced by a normalized integral 
lo^ 2 s ^ n Jo* over ^ ne hemisphere of axis orientations (by symmetry the other 
hemisphere yields an identical result) and a sum over the measurement outcomes m, 

Pi7 ] = 2J ^r r /2 Bini*W I*"** E PrnV,v)D*o(<p,#,0)1% m , (5) 
where the Stern-Gerlach probability distribution along a quantization axis ($, if) is 



given by the diagonal elements p mm of (2a I in the rotated frame, 

Pnfrtp) E [^(M.OIV^ ( 6 ) 
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Using the orthogonality relations of Clebsch-Gordan coefficients, 



E 



, j mm , j mm _ c /~\ 
l m Vo ~ d kk', {() 



and spherical harmonics, 

2 S in^^ 27r d^[^ (^^,0)]*^ (^^0) = (8) 

it is easy to show that indeed p^. — Pkq, proving the validity of the reconstruction 
method in the limit of infinitely many homogeneously distributed Stern-Gerlach 
experiments. 

In the more experimentally relevant case of a finite data set, the literature on 
the two-dimensional inverse Radon transform by filtered backprojection [29 indicates 
that excellent results can still be recovered, albeit with aliasing artifacts present to 
some degree. As a rough estimate, if Stern-Gerlach experiments are performed only 
along certain quantization axes spaced by an average angle Arj, then the reconstructed 
partial waves of the Wigner function become unreliable for k > fc max = ir / Ar/. Further, 
if the number M of measurements is much less than the number of degrees of freedom 
(fcmax + I) 2 , then the reconstructed coefficients P k q will be dominated by noise, in 
particular at large k. Both of these effects are mitigated in section [2~2] for the present 
reconstruction scheme. 

2.1. Accounting for fluctuations in the total angular momentum j 

We recall that for systems composed of many spin-1/2 components, such as two- 
component Bose-Einstein condensates, the total angular momentum j = (iV+ +JVi)/2 
often varies between single Stern-Gerlach measurements, as each such measurement 
requires the preparation of a new condensate. Instead of constructing a separate 
Wigner function for each occurring value of j, we notice that for k <C 2j the Clebsch- 
Gordan coefficients tu m depend smoothly on the total angular momentum j. This 
allows us to reconstruct the low-resolution part of the Wigner function even if j varies 
slightly between single Stern-Gerlach measurements. To this end we include the 
measured values of j in the data tuples, extending them to tp n , Cn, j n , m n ); the 
filtered backprojection formula is modified to 

M 

= (2fc + l)^c n ^ o (^,^,0)^ o m " m ". (9) 

n=l 

Again we refer to |Appcndix A| for a numerical method to evaluate this expression. 
The same smoothness of the Clebsch-Gordan coefficients at low k is used 



in section 2.2 to treat measurement uncertainties in both j n and m n in a 
perturbative manner in ^ . This is fundamentally different from a direct tomographic 
reconstruction of the Dicke matrix elements p mm > , where such uncertainties introduce 
large but correlated errors throughout the density matrix and make such a perturbative 
treatment impossible. 

2.2. Measurement uncertainties and high-k damping 

It is natural to assume that M uncorrelated experimental measurements can only serve 
to reconstruct M coefficients pk q , suggesting an upper limit fc max ~ VM (assuming 
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again a homogeneous distribution of quantization axis orientations) . For larger values 
of k the angular power spectrum |30| 

ft 



a 



(fb P ) 



2k 



(fbp)|2 
1 



(10) 



q=-k 



tends to acquire large fluctuations because of insufficient experimental data (see 
figure [3] for an example). However, simply cutting the reconstruction off at fc max 
is unsatisfactory because it disregards that some useful information is still present 
in these high-fc partial waves. A more natural cutoff is introduced through the k- 
dependent sensitivity to experimental uncertainties. Assuming experimental variances 
of (N$) - (A^) 2 = (Nf) - (Ni) 2 = a 2 N , we find that the uncertainties of (j 2 ) - (j) 2 = 
(m 2 ) — (m) 2 — <t%/2 (with no covariance, (jm) — (j)(m)) yield a leading order 
damping of the Clebsch-Gordan coefficients 

2 



\ l k0 I 



,'jmm 



' N 



-k(k+l) 



(11) 



2j(2i-l) 

The rotation matrix elements are damped similarly: if the pointing direction of the 
quantization axis Q = (■#, if) has an uncertainty of ern "C 1 {in terms of the expectation 
value of the angle 77^0' between the ideal axis orientation f2 and its true experimental 
value fi' we define CTq = (sin 2 rjaa. 1 ) = (1 — [cos -d cos$' + sin$sin$'cos((^ — if')] 2 )}, 
then for large k we find the rotation matrix elements to be damped as 



(£>*,(¥>, 0,0)) ™D k q0 (ip,#,0)exp 



4 



k(k + l) 



(12) 



If ct^v and (To are equal for all measurements, the linearity of ([£]) yields a simple 

smoothing p kq i-> Pkqe~ ak ^ k+1 ' with a — 2 j(2j-i) ~^~^4 L - ^ n tms wa ^' l nese tw0 damping 
formulas ( 11|12 ) cut off the reconstruction at large k in a natural and smooth way. 



2.3. Assembling the Wigner function 

Inserting the resulting coefficients (JoJ) into the form of the Wigner function ([lj we find 
the tomographically reconstructed Wigner function 

M 

W (ihv) (d,f) = 5^c, 
71=1 

M 

= c n Zj n ,m n [cos $ cos i? n + sin 1} sin $ n cos(y — <p n )]> (13) 
n=l 

where the contributions can be simplified to 

E jm (x) = -j= J2(2k + I) 3 / 2 t^ m P k (x). (14) 

v k=0 

As is to be expected in spherical symmetry, the contribution of an individual Stern- 
Gerlach measurement (see figure [T]) depends only on the relative angle cos T]cia n = 
cos $ cos i? n + sini? sini?,! cos((/9 — <^ n ) between the quantization axis orientation £l n = 
(i} n , <p n ) of the measurement and the point fi = (1?, if) on the Bloch sphere (figure [2]). 
Similarly to technical implementations of the planar inverse Radon transform |29j , the 
Wigner function is thus assembled from additive contributions due to the individual 



2j k 

EE 

fc=0 q=-k 



(2k + l)D k q0 (f n ,$ n , 0)Y kq (V, f)t 3 k l 
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Figure 1. Contributions S^ m (cos?7) to the Wigner function | |14| for j = 20 and 
m = — 20 . . . + 20. All curves have been divided by 100 and offset vertically by 
m. The bold curve for m = +16 is used in figure [2] Notice that the m = ±j 
contributions have lower spatial resolution (Arj ~ l/v(?) than those with m ~ 
(A.T] ~ 1/j); see section [3] 




Figure 2. Contribution to the Wigner function l |14| for j n = 20 and m n = 16 
(see figureJTJ; colors as in figure|4]but scaled to the maximum value of +163. The 
contribution H2o,i6( cos? ?) depends only on the angle rj between the quantization 
axis Q n and the direction Q in which the Wigner function is measured. 
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Stern-Gerlach measurements (see figure [2] for an example). The constructive or 
destructive interference of these contributions is what yields the reconstructed Wigner 
function (see figureB. The spatial resolutions of the Ej m (cosr]) ultimately determine 
the spatial resolution of the reconstructed Wigner function: if the Wigner function 
is composed predominantly of contributions with m n ~ ±j n its angular resolution is 
limited by that of a coherent state, Arj > if on the other hand the majority 

of contributions has m n ~ the resolution can be significantly higher, A77 > 
We make use of this observation in sections [3] and |4j where a spin-squeezed state is 
reconstructed and the increased spatial resolution is critically important. 



2-4- Positivity of the density matrix 

It is well known that only positive semi-definite density matrices represent valid 
quantum-mechanical states of a system pQ . Unfortunately, the filtered backprojection 
method ^ does not assure that the reconstructed p is positive semi-definite when used 
with a finite and noisy data set. For the purpose of displaying the Wigner function 
graphically, this is of no concern (see figure [4]) ; however, when the tomographically 
reconstructed coefficients p^^ are used in quantitative calculations (see section 41) 
positivity can be crucial. This is a similar problem as the requirement for a positive 
absorption density in medical computed tomography (CT) imaging |29j . It is also 
present in many quantum-state reconstruction schemes, and has been discussed 
extensively in the quantum tomography literature [T] . 

We do not offer a solution for assuring the positivity of the reconstructed density 
matrix. Here we merely point out that in other reconstruction schemes, such as 
maximum- likelihood estimates [3T], the ansatz p = T^T forces the density matrix p to 
be positive semi-definite; but a direct tomographic reconstruction of T similar to ^ 
is currently lacking. 



3. Quantization axes lying in a single plane 

When the spin-j system's quantum-mechanical state is fairly localized on the Bloch 
sphere, not every choice of quantization axis orientation has the same potential for 
extracting information about the state. When the axis is close to parallel to the 
state, most Stern-Gerlach measurements will yield \m\ j, with a limited angular 
resolution ~ 1/yfj given by the size of a coherent state on the Bloch sphere [26]. If 
the axis is close to perpendicular to the state, on the other hand, the distribution of 
measured values m represents the structure of the state's Wigner function much more 
accurately, with an angular resolution ~ 1/ j. This difference in scaling of the angular 
resolution, visible in figure [T] suggests that for large j it may be advantageous to 
focus on performing Stern-Gerlach measurements with quantization axes in a plane 
perpendicular to the quantum state, instead of covering the entire hemisphere of axis 
orientations. As a consequence much fewer measurements are needed, and we can get 
much more rapid convergence of the reconstruction in practice. But it is not a priori 
clear that this restriction of the quantization axes to a single plane has the potential 
for reconstructing the full quantum-mechanical state of the system. 

As it turns out, a modification to the "filter" function in ^ results in a full 
reconstruction of the mirror-symmetric part of the Wigner function. Defining the 
coordinate system such that the state is localized near the +z axis and all quantization 
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^ P) = (^)|(^)^E^o(^|,Q)4o m " m ": (15) 

71=1 

where (a) n = T(a + n)/T(a) is a Pochhammer symbol, and (a) 1 ~ ^fa — l/(8y / a). 

We again prove this reconstruction in the infinite-data limit. In the case of 
a homogeneous distribution of all azimuthal axis orientation angles ip we use the 
relationship 



5 gg > 



if k + q even 



(16) 

if k + q odd, 

which remains true in the experimentally more relevant case of a finite number A of 
equally-spaced axis orientations (replacing i Jq dtp i— >• -jj X^o* w ith f — an /A) as 
long as k < A. Together with ([6| and ^ we thus find that 




g -V 2 2 ^ J Q u ^ VmK^rl^qOWi ^' J k ° 

if fc + q even , 

(17) 

if k + q odd. 

Thus in the infinite-data limit such an in-plane reconstruction exactly determines the 
coefficients pk q for which k + q is even, while giving no information on the coefficients 
for which k + q is odd. Since the parity of k + q is the z •<->• — z reflection parity of 



the spherical harmonics Yfc g ($, tp), the in-plane formula ( 15 1 reconstructs the positive- 
parity component W + (d, (p) of the Wigner function W{$, ip) = W + {&, (p) + W~ ($, ip), 
with W ± (ir — fl,(p) = ±W ± (■&, ip). If we know from other measurements that the state 
is fully localized on the "northern" Bloch hemisphere (z > 0) , then the correct Wigner 
function is 

, , \2W + {d,ip) if Q < i? < f 

V Y> \ if f < ■& < IT, ^ ' 

which has the decomposition 

^p^n) = r sinMl? [** tyYMwwwM = J2 t L'pS p,p) (19) 

Jo Jo k>=0 

in terms of the overlap integrals T^ fe , given in Appendix B[ We conclude that the data 



acquired by Stern-Gerlach measurements with quantization axes lying solely within a 
plane are sufficient for an exact reconstruction of the Wigner function. 

3.1. Measurement uncertainties and high-k damping 



Measurement uncertainties can be introduced in ( 15 ) in the same way as in section [272] 
However, in an in-plane measurement series we can additionally separate out the 
azimuthal axis orientation uncertainty: since the rotation matrix elements D^ (ip, 0) 
are proportional to e~ lqv , a variance (ip 2 ) — (ip) 2 = u 2 ^ leads to a damping 

0)) = D* J ( VJ> | l 0)exp(-ig 2 4). (20) 
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Figure 3. Angular power spectrum \10\ of the reconstructed Wigner function 
of figure [1] Without damping (□) the power in modes k > 70 is too large and 
dominated by noise and aliasing effects; experimental uncertainties damp the 
angular power at large k in a natural way (•, see sections |2,2| and |3 . 1 \ . Odd-fc 
modes contain less power than even-fc modes because of the approximate point 
symmetry of the Wigner function (see figure Hb . 



4. Demonstration with experimental data 

In this section we reconstruct the Wigner function from a data set describing ensembles 
of N — 1250(45) atoms acquired in our group [BJ. In contrast to [BJ we rotate 
the coordinate system such that all quantization axes lie in the xy plane and the 
state is localized around the +z axis; in this way the procedure of section [3] can be 
employed directly. The data set consists of three experimental runs spanning different 
ranges of <p with different angular resolutions, owing to the fact that the need for 



homogeneity in ip for the filtered backprojection algorithm ( 15 ) was not known at the 
time of data acquisition. We use weights c n adjusted such that the weighted density 
of Stern-Gerlach measurements is as close to homogeneous as possible over the range 
ip = . . . 7r of azimuthal quantization axis orientations. As discussed in section [3] 
the planar arrangement of quantization axis orientations leads to a Wigner function 
which is peaked along both the +z and — z directions, featuring two identical copies of 
the quantum state. An additional Ramsey experiment [BJ was used to experimentally 
determine the correct location of the state on the northern (z > 0) Bloch hemisphere. 



High-fc damping (section 3.1) is achieved with an experimental uncertainty of 
(jjv ~ 11 atoms ( |11| and with an experimental error model dominated by phase noise: 
a v w <jp h sin(|^|)7v2 in (20), with phase noise amplitude cr p h = 8.2° [B]. In figure^ 



the effect of this damping is shown to be crucial for partial waves k > 70. 

The resulting reconstructed Wigner function is shown in figure [4] The high- 
frequency artifacts in the Wigner function far from the central peak are due to 
incomplete destructive interference of the contributions from the individual Stern- 



Gerlach measurements (see section 2.3). We expect that a more complete data set, 
including more quantization axis orientations, will lead to a smoother Wigner function 
at large angles 
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Figure 4. Reconstructed Wigner function using the data set described in 
section]!] The Wigner function takes values from —0.93 to +3.54. The coordinate 
system is rotated from [B] (see text). 



4-1- Spin-squeezing measurement 

We demonstrate the quantitative use of the reconstructed Wigner function by 
estimating the amount of spin squeezing in the system. Given a set of reconstructed 
Wigner function coefficients P k q we can calculate the probability distribution for the 
angular momentum projection quantum number onto any quantization axis orientation 
($,tp) from In principle the variance V(#, p) = (m 2 )($, <p) — \{m) (i?, tp)] 2 measures 
the amount of spin noise obtained experimentally. The expectation values of small 
integer powers of the projection quantum number m depend only on the low-fc 
components of the Wigner function, which are particularly insensitive to experimental 
noise flll|12|20| ); in particular, 



3 

E 



12 



J2 [Dla(v,#,0)]*Pi q 



m=-j 



. 7 (j + l)v / 27+T 



-POO 



(2j - 1) 5 



180 



]T [D 2 q0 (ip^,o)r P2q . 



(21) 



In figure [5] we plot the resulting variances for quantization axes in the xy plane, 
and compare them to a coherent state centered on the +z axis. In the presence of 



imaging noise the variance of such a coherent state is given by (21) with p k 



(con) 
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Figure 5. Normalized variance V = (m 2 ) — (m) 2 as a function of azimuthal 
quantization axis orientation ip. Open circles show variances calculated directly 
from Stern— Gerlach experiments along a given quantization axis [6]. The dashed 
line was calculated directly from the coefficients p^ P ' F ^ through (12 ll . The solid 
line shows the results of Gaussian fits (figure [6]l to the probability distributions 
Pm(J,v) given in tel. As in jjj we first subtract the experimental noise 
with <tjv = 11) fromthe calculated variances, and then divide by the variance of 
a coherent state, V co h = (j')/2 with (J) = 630 [see < |22[ |] , A spin-squeezed state is 
characterized by negative values (in dB). 
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Figure 6. Probability distributions p m of the projection quantum number m 
along the minimum- variance axis (ip B = 6.7°) and the maximum- variance axis 
(tp = -88.0°) of figure^ assuming j = 629 [found from p 00 = {(2j + l) -1 / 2 } 
0.02818]. Negative values of p m indicate that the reconstructed density matrix 
(as shown in figure [2| is not positive semi-definite and therefore does not strictly 
represent a physical state (see section |2.4| l. Gaussian fits used for figure [5] are 
shown as continuous lines. 
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(m2) (co h) = jV + 1) _ m-i) e -^ = i±^L +(D{(J%/n (22) 

3 6 2 

In figure [5] the experimental variance and the coherent-state variance are compared 
without imaging noise, i.e., the leading-order imaging noise contribution er^/2 is 
subtracted from the experimental variance before comparison with the variance of 
a noise-free coherent state. We notice that the resulting variance of the reconstructed 
state (dashed line in figure [5| is much larger than what was determined directly from 
the variances of the Stern-Gerlach data sets along the different quantization axes 
(open circles) . We believe that this is a result of the lack of positivity of the density 
matrix (see section 2.4 1, owing to the finite and noisy data set used for its tomographic 
reconstruction. In fact, the probability distributions in figure [6] clearly show negative 
values, which strictly speaking render the reconstructed density matrix unphysical. 

As an alternative extraction method wc calculate the probability distribution 
Pmi^: f) along a given quantization axis through ([6| and fit it with a Gaussian curve 
(see figure [6]); the variance of this fit then serves as an estimate of V(i9,ip). In such 
a fit the positivity of the p m is no longer a crucial ingredient. In figure [5] we show 
that this produces results that are very close to the variances calculated directly from 
Stcrn-Gcrlach experiments along the various quantization axes. The deviations close 
to the squeezing maximum (<^ s 6.7°) result from the fact that the reconstructed 
Wigner function contains contributions from all measurements, and therefore the 
extracted variance along a given quantization axis may be contaminated. Nonetheless 
the reconstructed Wigner function delivers a very concise picture of the structure 
of the multiparticle state, even for a data set with a non-uniform distribution over 
quantization axis orientations and with fluctuating values of j n . Further, with our 
method the variance V(i?, ip) can be calculated along any quantization axis orientation. 

In practice any proof of spin squeezing will not proceed through the reconstruction 
of the Wigner function followed by either a fit to the projection (JsJ) or a direct study 
of the projection noise (21). Instead, once the direction of squeezing tp a has been 
determined, a full Stern-Gerlach experiment will be performed along this axis in 
order to directly estimate the probability distribution p m (<£ s ), as in [B] and in figure [5] 
(circles). In this way, problems associated with the positivity of the reconstruction 
(section |2.4[ ) and with the influence of experimental data and noise from directions 
ip =/= (p s are strictly eliminated. 

In future experiments providing data for the present tomographic reconstruction 
method, wc plan to perform Stern-Gerlach measurements along many more 
quantization axes, but with as little as a single measurement per axis. Further, wc 
will pay attention to cover the entire range of quantization axes uniformly [either the 
entire equator for ( 15 ) or the entire sphere for Q] in each experimental run. In this 
way, we expect to need only minimal data preprocessing before reconstructing the 
Wigner function, and will be able to use the acquired data in the most efficient way 
by using equal weights c n = l/M for all data points. We also expect that for such an 
improved data set the variance of the simple estimate given by (21 ) will be closer to 
that of the quantum-mechanical state. 
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5. Conclusions 

We have presented a simple method for a tomographic reconstruction of the Wigner 
function of a spin-j system, applicable even to experimental settings where j is large 
and fluctuates between measurements. While the general procedure ([9| requires 
Stern-Gerlach type measurements spread uniformly over all possible quantization 



axis orientations, a more specialized and faster procedure ( 15 1 determines the Wigner 
function using only a single plane of quantization axis orientations. We have shown 
that this latter procedure is capable of reconstructing the Wigner function of a spin- 
squeezed state from a recently published experimental data set [5J. 
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Appendix A. Numerically evaluating Clebsch Gordan coefficients 



We have used a recursion relation to evaluate the Clebsch-Gordan coefficients 

-3,rn _ .jm-i 
'k — L k0 



r j, m _ iJmm _ (_!)i-m^ m; ^_ m | fc;0 ) from 3 : 



4j + l 

, 2 j - k 

2^'+ 1 /2 \\ (2j + l)_ 



T k ~ 



., j i _ k(k + 1) \ 



2j 



k 



, = 2j(j + l)-2(m+l) 2 -fc(fc + l) T . hm+1 _ j(j + i)-( m + i)( m + 2) rh m+2 
k j(j + 1) - m(m + 1) k j(j + 1) - m(m + 1) k 

ri'~ m = (-l) k 4' m - (A.l) 
This procedure is numerically stable even at very large values of j and k. 

Appendix B. Hemispherical overlap integrals of spherical harmonics 

The hemispherical overlap integrals of the spherical harmonics are |33) 

? q kk ,=2 f ' sintfdtf f d¥>y£O?,¥>)5W0,¥>) 
Jo Jo 
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= < 



f 1 if k = k' 

( 1 ) h ^2«- k - ± ^ 1 V / (2fc + l)(2fc / + l) 
1 ' (k-k')(k + k' + 1) 

I (fc - q)\(k' -q~)l (k' + q)ll(k + q- 1)!! 
X Y {k + q)\{k' + q)\ (^z|=l)!(fc=2)] 

if fc — g even and k' — q odd 



(-1)^— 2« 



v /(2fc + l)(2fc' + l) 



(B.l) 



(k - k'){k + k' + 1) 

. / (fc - g)!(fc' - g)[ (fc + q)\\{k' + g - 1)!! 
X Y {k + q)\{k> + q)\ (^l)l(^y. 



if fc' — g even and fc — g odd 



otherwise. 
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